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ABSTRACT 

Using observations of pulsars from the Parkes Pulsar Timing Array (PPTA) project 
we develop the first pulsar-based timescale that has a precision comparable to the 
uncertainties in international atomic timescales. Our ensemble of pulsars provides an 
Ensemble Pulsar Scale (EPS) analogous to the free atomic timescale Echelle Atom- 
ique Libre (EAL). The EPS can be used to detect fluctuations in atomic timescales 
and therefore can lead to a new realisation of Terrestrial Time, TT(PPTAll). We 
successfully follow features known to affect the frequency of the International Atomic 
Timescale (TAI) and we find marginally significant differences between TT(PPTAll) 
and TT(BIPMll). We discuss the various phenomena that lead to a correlated signal 
in the pulsar timing residuals and therefore limit the stability of the pulsar timescale. 

Key words: pulsars: general — time 



1 INTRODUCTION 

Atomic frequency standards and clocks are now the basis 
of terrestrial time keeping. Many countries distribute a lo- 
cal atomic timescale. These are combined by the Bureau 
International des Poids et Mesures (BIPM) to form Inter- 
national Atomic Time (or Temps Atomique International, 
TAI) which is published in the form of differences from the 
national timescale^]. TAI is the basis for both Coordinated 
Universal Time (UTC), used for the dissemination of time 
signals, and Terrestrial Time (TT). TT is formed by ref- 
erencing individual clocks to the Earth's geoid. Through- 



The differences between TAI and various other timescales 
can be obtained from the "Circular T" publication avail- 
able from |http : / /www . bipm . org/en/scientif ic /tai/ 1 The dif- 
ference between TT(BIPM) and TT(TAI) is provided at 
|f tp : //tai . bipm ■ org" /TFG/TT (BI PM) | 



out this paper, we refer to TT(TAI) as terrestrial time re- 
alised by TAI. Once published, TAI itself is never revised, 
but the BIPM publishes another realization of TT which is 
computed every year and labelled TT(BIPMYY), where YY 
corresponds to the year of the most recent data used. For 
instance, in this paper we refer to TT(BIPMll) as the most 
recent post-corrected realisation. 

The difference between TT(BIPMll) and TT(TAI) is 
shown in the top panel of Figure [1] and clearly shows a drift 
between the time standards of ~ 5 /j.s since 1994. The stabil- 
ity of TAI is obtained from a large number of atomic clocks 
whereas the accuracy of TAI is set from a few primary fre- 
quency standards (Arias, Panfilo & Petit 2011). Initially, 
the free atomic timescale Echelle Atomique Libre (EAL) 
is produced from the weighted average of the timescales 
of several hundred atomic clocks around the world. This 
timescale is not in accord with the second as defined in the 
International System of Units (SI). Therefore, to form TAI 
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Figure 1. The top panel shows the difference between 
TT(BIPMll) and TT(TAI) since the year 1994. The bottom 
panel shows the same, but after a quadratic polynomial has been 
fitted and removed. 

(which does conform to the SI second) from EAL, various 
frequency adjustments are necessary. These are determined 
using primary frequency standards. Frequency adjustments 
are generally made slowly, a process referred to as "steer- 
ing". In 1996, a decision was made to change the realization 
of the SI second that resulted in a frequency shift of about 
2 x 10 ?14 . That shift was progressively introduced into TAI 
over a period of two years. As TAI itself is never retroac- 
tively corrected, only the post-corrected versions of TT, e.g., 
TT(BIPMll), have the earlier data corrected. This leads to 
the "bump" that we observe in Figure 1 around the year 
1998. 

Although numerous clocks are used in forming TAI and 
there is continuous development of atomic clocks, stability 
over decades is difficult to measure and maintain. It is there- 
fore desirable to have an independent precise timescale valid 
on such long intervals. In this paper, we describe the devel- 
opment of such a timescale based on the rotation of pulsars. 

Radio pulsars are rotating, magnetised neutron stars 
that radiate beams of electromagnetic waves. For a fortu- 
itous line of sight to the pulsar, these can be observed at the 
Earth as pulses. The pulse times of arrival (ToAs) from the 
brightest and fastest-spinning pulsars can be measured with 
a precision of ~ 100 ns in an observation time of ~ 1 hour. 
This precision is significantly worse than that obtainable 
from atomic clocks, but, in contrast to individual clocks, 
can be maintained for a very long time. We note that a 
pulsar-based timescale provides: 

• an independent check on terrestrial timescales using a 
system that is not terrestrial in origin. 

• a timescale based on macroscopic objects of stellar mass 
instead of being based on atomic clocks that are based on 
quantum processes. 

• a timescale that is continuous and will remain valid far 
longer than any clock we can construct. 

In order to develop a pulsar-based timescale, all phe- 
nomena affecting the pulse ToAs must be taken into ac- 
count. These are incorporated into a "pulsar timing model" 
that contains the pulsar's astrometric, rotational and or- 
bital parameters, the effects of the interstellar medium and 



the motion of the Earth about the solar system barycen- 
tre. Timing residuals are the difference between the arrival 
times converted to the solar system barycentre and predic- 
tions of those times based upon the timing model (see e.g., 
Edwards, Hobbs & Manchester 2006 for details). Non-zero 
residuals can result from an incorrect conversion from the 
measured ToAs to barycentric arrival times. Our ability to 
convert to barycentric arrival times relies, for instance, upon 
the accuracy of the solar system ephemeris. Many pulsars 
also display irregularities in rotation and changes in pulse 
shape that make timing difficult (e.g., Lyne et al. 2011 and 
references therein). A subset of pulsars, the "millisecond pul- 
sars", have shorter pulse periods and much more stable ro- 
tation than the "normal pulsars". However, precise observa- 
tions of millisecond pulsars show some unexplained timing 
irregularities which we refer to as "timing noise" . 

Some of the variations in the timing residuals are caused 
by processes that are correlated between different pulsars. 
These can be identified by observing an ensemble of pul- 
sars, a so-called "Pulsar Timing Array" (PTA, e.g., Foster 
& Backer 1990). Errors in the terrestrial time standard will 
introduce exactly the same signal in the residuals for each 
pulsar. In contrast, errors in the planetary ephemeris used in 
the timing analysis will induce timing residuals which have a 
dipolar signature on the sky and gravitational waves prop- 
agating past the pulsar and the Earth will induce timing 
residuals with a quadrupolar signature. As shown later in 
this paper, it is not possible to obtain an unbiased estimate 
of the time standard errors simply by forming a weighted 
average of the timing residuals for different pulsars. This is 
because of the coupling between the timing model for each 
pulsar and the measurement of the correlated signal as well 
as the differing data spans for each pulsar. 

In this paper we analyse data from the Parkes Pulsar 
Timing Array (PPTA) project (Manchester et al. 2012) to 
develop a pulsar-based timescale which we label an Ensem- 
ble Pulsar Scale (EPS). This scale has similarities to the 
free atomic timescale EAL. The frequency of EAL needs to 
be steered using primary frequency standards to realise a 
timescale based on the SI second. Similarly, since the intrin- 
sic pulsar pulse periods and their time derivatives are un- 
known for the pulsars in a PTA, the EPS is not an absolute 
timescale and it must be "steered" to a reference timescale 
which conforms to the SI. This is achieved by first forming 
timing residuals for each pulsar with respect to the reference 
timescale, TT(TAI) in our case, and subsequently fitting a 
quadratic polynomial to the residuals. Fluctuations in the 
reference timescale with respect to the EPS can be identified 
and used to provide a set of corrections to that realisation 
of TT, thereby realising a new pulsar-based timescale. We 
refer to the timescale derived in this paper as TT(PPTAll). 
The bottom panel of Figure [T] shows the difference between 
TT(BIPMll) and TT(TAI) after a quadratic polynomial 
has been fitted and removed. It is this signal that we expect 
to see in comparing TT(PPTAll) with TT(TAI). 

Earlier attempts to develop a pulsar timescale have been 
made by Guinot & Petit (1991), Petit & Tavella (1996), 
Rodin (2008) and Rodin & Chen (20110. We will show be- 
low that, in contrast to our method, these earlier attempts 



2 Note that some authors (e.g., Petit & Tavella 1996; Rodin, 
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did not account fully for the effects of fitting a pulsar tim- 
ing model. They also have not been applied to high precision 
observations for a large number of pulsars. 

In §2 we describe the signal that is potentially measur- 
able using pulsar observations. §3 describes the observations 
used in this paper. §4 contains details of the method ap- 
plied. §5 presents the application of our method to actual 
data and contains a discussion on the result. §6 summarises 
the results. The algorithm presented here has been included 
in the tempo2 pulsar-timing software package (Hobbs, Ed- 
wards & Manchester 2006). Usage instructions are given in 
Appendix A. 



2 THE CORRELATED SIGNAL 

Pulsar timing models are based on the proper time, t psr , 
measured at the centre of the pulsar assuming that its grav- 
itational field is not present. Note that the actual time of 
emission of a pulse and its time of arrival at the solar sys- 
tem barycentre differ by the light travel time from the pul- 
sar, which, for this work, is assumed to be constant- The 
time of emission of a pulse from the pulsar, t^ SI , is therefore 
related to the observed ToA, t° bs , as 

C r = t° bs + A clk + A pc + A nc . (1) 

A c ik includes all the steps required to convert the measured 
ToA to barycentric coordinate time (TCB). The steps (listed 
below) for this correction are identical for different pulsars. 
Any error in this correction will therefore lead to timing 
residuals for different pulsars that are exactly correlated 
(i.e., Acik(t) will be identical for all pulsars). A pc represents 
steps in the processing that lead to timing residuals that are 
partially correlated between different pulsars. For instance, 
the correlation coefficient may depend upon the angle be- 
tween the pulsars. A nc represent corrections that are specific 
to a given pulsar and are not correlated between different 
pulsars. In addition to pulsar dependent effects these cor- 
rections include the effects of the interstellar plasma and 
radiometer noise. 

Errors in the timing system (A c ik) lead to timing resid- 
uals that are correlated between multiple pulsars. A c ik can 
be separated into various components as 

A c ik = Atrp + Att + Atcb- (2) 

Atrp represents the time delay between the topocentric refer- 
ence point of the telescope and the time tagging on the out- 
put data. Most of these delays are constant (for instance, at 
the Parkes Observatory there is an approximate time delay 
of 600 ns from the receiver to the backend instrumentation 
that records the signal). However, each backend instrument 
also has an effective delay and such delays differ for each 
instrument. At the Parkes Observatory these delays can be 
tens of microseconds. A method for measuring and correct- 
ing for such delays is described in Manchester et al. (2012). 
It is not currently clear how stable these delays are. Initial 

Kopeikin & Ilyasov 1997) have considered using the orbital pa- 
rameters of binary pulsars to provide a pulsar-based timescale. 

3 We note that all pulsars have a radial velocity. The effect of this 
velocity is to change the observed pulse frequency by an effectively 
fixed amount. 



studies have suggested that, for some of the backend instru- 
mentation, variations at the 10-100 ns level may be occur- 
ring between observing sessions at the Parkes Observatory. 
However, inaccuracies in measuring these delays generally 
leads to step-changes when a new observing system is com- 
missioned, or adds high-frequency noise if the instrumental 
time delays randomly change between observations. These 
effects are therefore different to the secular drifts expected 
from errors in terrestrial timescales (see Figure [T]). 

Att is the time difference between the observatory 
clock and a reference implementation of TT such as 
TT(TAI). In order to determine the approximate uncer- 
tainty in Att, we compare two independent techniques. The 
first method, which is used for our standard data processing, 
converts from the Observatory time standard to the Global 
Positioning System (GPS) time standard and uses tabulated 
corrections from the GPS system to terrestrial time. The 
second method uses a GPS Common View system to trans- 
fer the Observatory time to the Australian time standard, 
UTC(AUS). Tabulated corrections are subsequently used to 
convert from UTC(AUS) to terrestrial time. We determined 
the difference between these two techniques every 10 d dur- 
ing the year 2011. The rms difference between the two meth- 
ods is 8.8 ns implying that the precision of the time transfer 
is of this order. 

To obtain barycentric arrival times we have to convert 
from TT to TCB. Conversion from TT to TCB is carried 
out using a time ephemeris described by Irwin & Fukushima 
(1999). In their paper, it is shown that this time ephemeris 
is known to better than 5 ns and therefore any errors will 
not significantly affect the pulsar timing residuals. 

A pc in Equation [1] represents corrections that, if they 
are not known with sufficient precision and accuracy, can 
lead to timing residuals for different pulsars that are par- 
tially correlated. All ToAs are corrected for the geometrical 
time delay between the Observatory and the solar system 
barycentre. This is carried out using the Jet Propulsion Lab- 
oratory DE421 solar system ephemeris (Folkner et al. 2008). 
Any inaccuracies in this ephemeris will lead to timing residu- 
als whose amplitude depends upon the position of the pulsar 
with respect to the ecliptic plane. For two pulsars that are 
close together on the sky, ephemeris errors will lead to tim- 
ing residuals that are correlated. However, for widely sep- 
arated pulsars, ephemeris errors will lead to anticorrelated 
residuals. At present, errors in the mass of Jupiter and Sat- 
urn are the most likely observable effects (Champion et al. 
2010); it is very likely that pulsar observations will improve 
our knowledge of these masses in the next decade. We will 
discuss ephemeris errors in more detail in Section [5] 

The main scientific driver for pulsar timing array ob- 
servations is the possibility of detecting gravitational waves. 
These result in variations in the timing residuals with a 
quadrupolar signature. The phenomenon thought to be the 
most likely to be detected is an isotropic, stochastic, grav- 
itational wave background (Hellings & Downs 1983). Such 
a background will induce a correlation of —0.15 < £ < 0.5 
between pulsar pairs depending upon the angle between the 
pulsars. For our sample of pulsars, the mean |£| is 0.15, mean 
C is 0.02 and the maximum £ is 0.42 for PSRs J1730-2304 
and J1744— 1134. We discuss the possibility that a gravita- 
tional wave signal could be misidentified as a clock error in 
Section [5] 
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Figure 2. The timing residuals used in this analysis referred to 
TT(TAI). The label on the right-hand side gives the pulsar name 
and the total range of the residuals for that pulsar. 



3 OBSERVATIONS 

The data used here are the extended Parkes Pulsar Timing 
Array (PPTA) data set that is described in Appendix A of 
Manchester et al. (2012). This data set included observa- 
tions taken as part of the PPTA project that commenced 
in the year 2005 along with earlier observations published 
by Verbiest et al. (2008, 2009). All observations were ob- 
tained using the Parkes 64-m radio telescope. Typical ob- 
servation durations were 1 hr. The observing system has im- 
proved significantly since the earliest observations and so the 
ToA uncertainties have generally decreased with time. We 
have measured the majority of the timing offsets between 
the different observing systems. Hence, we have been able 
to remove most, but not all, of the arbitrary offsets from 
the timing model that were included in the Verbiest et al. 
(2008, 2009) analyses (details are given in Manchester et al., 
2012). As the timing residuals for PSR J1939+2134 exhibit 
timing noise that is currently uncorrectable and is at a level 
significantly higher than the white noise level, we have not 
included this pulsar in our analysis. 

Observations since 2005 have been corrected for vari- 
ations in the ionised interstellar medium using multi- 
frequency observations, but observations prior to 2005 could 
not be corrected because adequate multi-frequency observa- 
tions were not made. Dispersion measure fluctuations are a 



significant noise source for many of the pulsars in the PPTA 
(You et al. 2007), so our inability to correct them in the 
early observations leads to larger error bars on the estimated 
TT(PPTAll) before 2005. 

Timing residuals were formed using the tempo2 soft- 
ware using the JPL DE421 solar system ephemeris and re- 
ferred to TT(TAlj3 The timing residuals for the 19 pulsars 
are shown in Figure [2] and a summary of our data sets is 
provided in Table [T] where, in column order, we provide: 1) 
the pulsar name, 2) pulse period, 3) dispersion measure, 4) 
weighted rms residual, 5) unweighted rms residual, 6) me- 
dian ToA uncertainty, 7) data span, 8) number of observa- 
tions, 9) date of first observation and 10) date of most recent 
observation. We emphasise that the following properties of 
the data set must be accounted for in the analysis: 

• Each pulsar has a different data span. For some pulsars 
data exist from 1995 onwards. For other pulsars only 6 — 8yr 
of data exist. 

• The data sampling is irregular with the more recent 
data being more uniform than earlier data. 

• Very few observations were made around the year 2000. 
Only PSR J0437— 4715 provides significant data around this 
time. 

• The ToA uncertainties are variable. They generally de- 
crease with time as new instruments were commissioned. 
However, pulsar scintillation also leads to significant varia- 
tions in the uncertainties. It is also common that the uncer- 
tainties underestimate the white noise present in the data. 
We account for this by including scaling factors that increase 
the error bars. In Table[2]the median ToA uncertainty is de- 
termined without these extra scaling factors. 

• Timing noise is observed in many of the data sets. 

• The rms timing residuals vary widely. The small- 
est weighted rms residual is 0.23 fis for PSR J0437-4715, 
whereas the largest is 5.1 /xs for PSR J1045— 4509. Because 
of red noise in many of the data sets, these rms values are 
often larger than the typical ToA uncertainties. 



4 METHOD 

To include clock errors in the timing model we need a func- 
tion that describes the clock error, A c (t), at any time, t, 
during the data span, in terms of a small number of param- 
eters. We want to avoid imposing any more structure on this 
model than necessary. We tried two approaches: a Fourier 
series and a set of equally spaced samples with an interpo- 
lation mechanism. Both provided adequate results, but we 
found the set of equally spaced samples provided better error 
estimates and more flexibility with the model parameters. 
We use linear interpolation between the samples which have 
spacing T s ; this is equivalent to a low-pass filter with a band- 
width of }lp = l/2T a . Regardless of the functional form of 
the model, the parameterisation of the clock error will have 
some covariance with the other parameters in the timing 
model. We implemented constraints on the clock parame- 
ters to minimise this covariance. These constraints zero the 



4 This is in contrast to the same data set described by Manch- 
ester et al. (2012). In that paper the data were referred to 
TT(BIPMll). 
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Table 1. Parameters for the pulsar timing residuals referred to TT(TAI). 
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Figure 3. Simulated timing residuals for three pulsars in the 
presence of the deliberate steering of TAI. The left-hand panel 
shows the pre-fit residuals for the three pulsars. The right-hand 
panel (which has a different y-scaling) shows the post-fit residuals. 



offset, linear and quadratic terms in A c (t), and those terms 
that represent a position error, parallax or proper motion. 
The changes to the least-squares-fitting algorithm that en- 
able these constraints are described in detail by Keith et 
al. (2012). As these terms are removed from the individual 
pulsar residuals they must not exist in the clock error be- 
cause they would be unconstrained. The modifications to 
the standard tempo2 least-squares-fitting procedure to al- 
low the A c (t) values to be fit globally to all these data sets 
were originally described by Champion et al. (2010). 

It would be simpler to form the weighted average of the 
timing residuals in order to determine the correlated signal. 
This is not useful because the resulting clock errors in the 
timing residuals for a particular pulsar will be modified by 



2000 



2005 



2010 



Year 



Figure 4. In all panels the dotted line represents TT(BIPMll)- 
TT(TAI) with a quadratic polynomial fitted and removed. The 
solid line in the top panel shows the result from a simple weighted 
average of the simulated timing residuals that are shown in Fig- 
ure [3] The data points in the central and lower panels are the re- 
sult of the algorithms described in this paper. In the lower panel 
the error bars account for the possibility of timing noise in the 
timing residuals for each pulsar. 



the fitting process that has been carried out for that spe- 
cific pulsar. In order to illustrate this effect, we simulate 
the timing residuals for three synthetic pulsars. Each pulsar 
has the same ToA uncertainty (50 ns) and is sampled every 
14 days. The first pulsar has continuous observations from 
the year 1994 to 2011. The second pulsar only has obser- 
vations until 2001 and the third pulsar has a gap of a few 
years around the year 2003. The observations are simulated 
assuming that TT(BIPMll) is perfect, but residuals are 
formed using the TT(TAI) timescale. These resulting pre-fit 
residuals therefore exhibit the differences between TT(TAI) 



6 G. Hobbs et al. 



and TT(BIPMll) along with 50 ns of additional white noise. 
A standard pulsar-timing fit is carried out for each pulsar 
(i.e., the pulse frequency and its first time derivative are fit- 
ted for). For the third pulsar, we also include an arbitrary 
phase jump between the early and late observations. The 
pre- and post-fit timing residuals are shown in left and right 
panels of Figure [3] respectively. The pre-fit residuals, in the 
left-hand panel, clearly are correlated and take the expected 
form shown in Figure [T] (note that due to the nature of the 
simulations these residuals are the inverse of the clock error 
shown in Figure [TJ. However, the fitting procedure signifi- 
cantly modifies the shape of the residuals and the post-fit 
residuals have a correlation coefficient significantly less than 
1. 

The average of our simulations is shown in the top panel 
of Figure 13] The average was simply calculated as the mean 
residual for each samplcfl Even though the average does 
show some of the features of the clock error, it does not 
model the clock errors perfectly because of the fitting of the 
pulsar timing models. For instance, the weighted average 
also leads to a step-change around the year 2001; this occurs 
when there is a change in the number of pulsars contributing 
to the average. Such a procedure will therefore only work if 
all the pulsars have an equal data span and the same timing 
model fits are applied to all pulsars. 

A c (t) values obtained using our method and their un- 
certainties are shown in the central panel of Figure[4] Clearly 
this procedure successfully models the features resulting 
from the steering of TAI. We note that the uncertainties on 
A c (i) remain relatively small between the years 2001 and 
2004 even though we have simulated data for only one pul- 
sar during this time. However, in reality, it is not possible 
to distinguish between clock errors and pulsar timing noise 
using data from a single pulsar. To account for both timing 
noise and clock errors we first assume that all the noise in 
a given pulsar's timing residuals is timing noise. We model 
the spectrum of this noise (using the SPECTRAlModel plu- 
gin to tempo2) and use a generalised least-squares-fitting 
procedure (Coles et al. 2011) to account for the timing noise 
when fitting the pulsar timing model and the A c (t) param- 
eters. The bottom panel in Figure U demonstrates how the 
error bars on A c (£) significantly increase when the noise for 
each pulsar is modelled as timing noise. If, as in this case, a 
significant clock error is measured, then these errors can be 
included in the timing procedure and the process iterated 
to determine true spectral models of the timing noise and 
hence the true error bars on A c (t). 

We emphasise that the A e (t) measurements and their 
uncertainties are not necessarily independent. The effect of 
fitting, irregular sampling, differing data spans and the lin- 
ear interpolation between adjacent grid points will all lead 
to correlated A c (t) values. The amount of correlation can 
be determined from the covariance matrix of the fit. 

In order to test our algorithm with realistic data, we 
formed simulated data with the exact sampling and ToA 



5 Note that previous work based on the average timing residuals, 
such as Rodin & Chen (2011), use a Wiener filter when averag- 
ing their data sets. We are able to obtain the mean residual for 
each sample because our simulations have equal weighting and 
identical sampling. 
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Figure 5. Results obtained for all 19 PPTA pulsars, but with 
simulated arrival times. In (a) only white noise and the clock 
error is simulated, (b) only white noise and c) only uncorrelated 
red noise. In (d) the same simulations are carried out as in (c), 
but PSR J0437-4715 is not included. 

uncertainties as in the observations of the 19 actual pulsars. 
Initially we added no timing noise, but as before, simulated 
the data using TT(BIPMll) and formed timing residuals us- 
ing TT(TAI). In the top panel of Figure [5] we demonstrate 
that our algorithm correctly recovers the expected signal. 
This demonstrates that the data sets and ToA uncertain- 
ties are such that our algorithm should correctly recover the 
irregularities in TT(TAI). 

We wish to confirm that our method does not incor- 
rectly lead to an error in TT if none exists. In Figure [5p 
we show the results after forming the timing residuals using 
TT(BIPMll). In this case, no clock error exists in the data 
(the timing residuals are purely white noise). We correctly 
find no significant A c (t) values. For Figure[SJ; we have simu- 
lated white data and added uncorrelated red noise to repre- 
sent timing noise. The resulting clock function does not show 
an unexpected large signal, but does show some correlated 
structure in the resulting data points. As discussed in more 
detail below, this is mainly caused by the timing residuals for 
PSR J0437-4715 dominating the fit. In Figure [SJI we repro- 
duce the analysis, but do not include PSR J0437-4715. The 
addition of timing noise therefore increases the error bars on 
A c (t), but does not lead to an incorrect measurement of a 
clock error. 



5 RESULTS AND DISCUSSION 

The top panel in Figure [S] indicates the sampling for each of 
the 19 pulsars. The solid line in the lower panel is the dif- 
ference between TT(BIPMll) and TT(TAI). This has been 
obtained by 1) sampling the expected signal at the same 
times as our measured clock function, 2) fitting a quadratic 
polynomial using an unweighted least-squares-fit and 3) re- 
moving this quadratic polynomial from the expected sig- 
nal. The data points in the Figure represent A c (t), the dif- 
ference between the pulsar-based timescale TT(PPTAll) 
and TT(TAI). All the recent data are consistent within 
2a with the expectation from TT(BIPMU), however some 
marginally significant differences are seen in the earlier data. 
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Figure 6. The top panel shows the sampling for the 19 pul- 
sars in our sample. The lower panel shows the difference between 
TT(BIPMll) and TT(TAI) as the solid line. The data points 
indicate the difference between TT(PPTAll) and TT(TAI). 
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Figure 7. The correlated signal caused by (a) errors in the Solar 
System ephemeris and (b) from one realisation of a gravitational 
wave background with dimensionless amplitude of 10 — 15 . The 
solid line indicates the expected correlated signal caused by the 
steering of TAI. 



A careful statistical analysis is non-trivial as 1) the error bars 
on the data points are correlated and 2) the clock errors are 
constrained so that they do not include a quadratic poly- 
nomial. The reduced-x 2 value obtained by comparing the 
expected clock signal with our data is 2.7. This simple sta- 
tistical test assumes that each data point is independent, but 
the value does indicate that there are no large discrepancies 
between the expected clock errors and the measurements. 

The most obvious discrepancies between our values and 
the expectation occur between the years 1995 and 2003. 
However, 1) our pulsar data set has sparse sampling around 
this time and has not been corrected for dispersion mea- 
sure variations and 2) the observed discrepancies would re- 
quire an error in the frequency of TT(BIPMll) of ~ 10~ 14 
whereas the uncertainty on this frequency is thought to be 
~lx 10~ 15 around the year 2003 (Petit 2003). This sug- 
gests that the discrepancies result from the determination 
of A c (t). There may be sufficient archival observations from 
other observatories to improve the clock error estimates dur- 
ing this period and thus to confirm or deny these possible 
errors in TT(BIPMll) and its estimated uncertainty. 

It is possible that errors in the solar system ephemeris 
could lead to correlated signals in the timing residuals. To 
see the maximum size of any such signal in our data, we 
have simulated observations using the same sampling and 
ToA uncertainties as the real data using the JPL DE421 
solar system ephemeris, but without any clock errors. We 
then processed the data using the earlier JPL DE414 solar 
system ephemeris. The resulting estimate of the "clock er- 
rors" are shown in the top panel of Figure [JJ The maximum 
deviation for recent data is < 100 ns. As we use the most re- 
cent ephemeris, DE421, for our analysis it is likely that the 
actual correlated signal caused by the planetary ephemeris 
is significantly smaller than this. 

In order to test whether a gravitational-wave back- 
ground signal could be mis-identified as an error in the ter- 
restrial time standard, we have simulated multiple realisa- 
tions of a gravitational-wave background (Hobbs et al. 2009) 
with a dimensionless strain amplitude of 10 -15 . This ampli- 
tude is typical of that expected for a background created 



by coalescing supermassive black-hole binary systems (e.g., 
Sesana, Vecchio & Colacino 2008). For each simulation we 
use the real sampling and ToA uncertainties as in the actual 
observations. The results from our algorithm for one realisa- 
tion are shown in the bottom panel of Figure [JJ This shows 
that, for current data spans, it is unlikely that such a signal 
will significantly affect the stability of the pulsar timescale. 
However, with increasing data lengths and with improve- 
ments in the ToA precision achievable, the gravitational- 
wave background could become a significant factor. In this 
case the clock estimation algorithm would need to be mod- 
ified to make it orthogonal to the gravitational wave back- 
ground. 

From our data, we therefore conclude that 

• the difference between TT(TAI) and TT(BIPMll) can 
be detected using pulsar data and that this difference, as 
expected, results from the deliberate steering of TAI. 

• there are no large unexpected errors in TT(BIPMll) 
over our data span. 

• the variations in TT(TAI) are at a significant level com- 
pared with the precision of current pulsar timing array ob- 
servations. We note that Guinot (1988) and subsequent pa- 
pers from the clock community have already pointed out 
that TT(TAI) is not suitable for high time precision pulsar 
experiments and that TT(BIPM) should always be used. We 
confirm that TT(BIPMll) is adequate for current millisec- 
ond pulsar timing experiments. 

Our results do not show the steering of TAI as clearly 
as in Figure [5^. This is mainly because the timing resid- 
uals for PSR J0437— 4715 dominate the data set: it is the 
only pulsar that was observed around the year 2000 and 
has a large number of observations and very small ToA un- 
certainties. However, the statistical properties of the tim- 
ing residuals for this pulsar suddenly change around the 
year 2006. Prior to this date, observations were made in 
the 20 cm band and have not been corrected for dispersion 
measure variations. After this date, the observations were 
made with new instrumentation, in the 10 cm band and the 
dispersion measure variations have been measured and re- 
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moved. Unfortunately, we do not currently have any method 
that can model highly non-stationary noise. The statistical 
model of the timing noise is applicable over the entire data 
set and, as such, is not optimal for any individual section. 
Future algorithmic developments may allow the effects of 
non-stationary noise to be included in our standard analy- 
sis procedure. Prior to 1996, our data sets are dominated by 
PSRs J1713+0747 and J1744-1134 which have significantly 
lower rms timing residuals than the other pulsars observed 
during that time. However, the effects of observations of just 
a few pulsars dominating the fitting procedures and the ef- 
fects of non- stationary noise will, in subsequent work, be 
mitigated by including observations of more pulsars from 
other observatories. 

The International Pulsar Timing Array (IPTA) project 
is a collaboration between three individual projects based in 
Europe, Australia and North America (Hobbs et al. 2010). 
The total number of pulsars being observed changes as new 
pulsars are discovered, but approximately 40 pulsars arc 
currently being observed. Of these, the ToA timing pre- 
cisions for ~30 should be better than 1/is. Production of 
high-quality data sets is ongoing and will soon lead to a sig- 
nificantly improved pulsar timescale. In the longer term, the 
Square Kilometre Array telescope (e.g., Cordes et al. 2004) 
should be able to observe many hundreds of pulsars with 
a timing precision of 100 ns or better. If, as expected, pul- 
sars are stable over long timescales at this level then such 
data sets should provide a long-term time standard that is 
competitive with the world's best terrestrial time standards. 



6 CONCLUSION 

We have developed a new algorithm for determining the 
correlated signal in the timing residuals for multiple pul- 
sars. Any errors in the reference timescale will lead to 
such a correlated signal. By comparing our measurements 
of pulse arrival times to TT(TAI) we have confirmed that 
we can recover the effects of the deliberate steering of 
TAI. We have not identified any significant discrepancies 
with TT(BIPMI I ), but have noted a marginal discrep- 
ancy between 1995 and 2003. Other phenomena, such as an 
isotropic, stochastic, gravitational wave background will also 
lead to a correlated signal in pulsar data sets, but we show 
that such phenomena are not likely to affect our results. 

In the future it is likely that pulsar data sets will be 
processed in a manner that will simultaneously identify: ir- 
regularities in the time standard; errors in the solar system 
ephemeris; and gravitational waves. By combining observa- 
tions from numerous telescopes, such future data sets will 
significantly improve on the pulsar-based timescale that is 
presented here. 
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APPENDIX A: SOFTWARE INSTRUCTIONS 

The algorithm described in this paper has been included in 
the CLOCK plugin to tempo2. Assuming that the user has a 
set of parameter and arrival time files (.par and .tim), the 
following procedure can be followed: 

• Obtain a spectral model for the timing noise in each 
pulsar and the corresponding covariance function: 

> tempo2 -gr spectralModel -f psrl.par psrl.tim 

(repeat for each pulsar). 

• Create a global parameter file (global . par) that con- 
tains the required realisation of Terrestrial Time and a set 
of IFUNC parameters that define the U grid points described 
in the text: 
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# global. par 
CLK TT(TAI) 
EPHEM DE421 
SIFUNC 2 2 
IFUNC1 49300 
IFUNC2 49600 

IFUNC22 55600 

where the IFUNC values range from before the earliest ob- 
servation to after the latest observation. Note that the 
SIFUNC 2 2 selects linear interpolation between the grid 
points (the first '2' on this line) and states that this pa- 
rameter should be fitted globally between the pulsars (the 
second '2' on this line). 

• The global fit must be constrained not to include an 
offset, linear or quadratic component. The following should 
be included in the parameter file for the first pulsar: 

CONSTRAIN IFUNC 

• The clock plugin can now be run. Typical usage is: 

> tempo2 -gr clock -fitfunc globalDCM -global 
global. par -f psrl.par psrl.tim -f psr2.par 
psr2.tim -f psr3.par psr3.tim ... 



